Thermodynamically-guided machine learning modelling for predicting the glass-forming ability of bulk metallic glasses

Glass-forming ability (GFA) of bulk metallic glasses (BMGs) is a determinant parameter which has been significantly studied. GFA improvements could be achieved through trial-and-error experiments, as a tedious work, or by using developed predicting tools. Machine-Learning (ML) has been used as a promising method to predict the properties of BMGs by removing the barriers in the way of its alloy design. This article aims to develop a ML-based method for predicting the maximum critical diameter (Dmax) of BMGs as a factor of their glass-forming ability. The main result is that the random forest method can be used as a sustainable model (R2 = 92%) for predicting glass-forming ability. Also, adding characteristic temperatures to the model will increase the accuracy and efficiency of the developed model. Comparing the measured and predicted values of Dmax for a set of newly developed BMGs indicated that the model is reliable and can be truly used for predicting the GFA of BMGs.

T g +T l ) could guide scientists to compare the glass-forming ability of alloy systems. Tripathi et al. 16 , combining the thermodynamics and principles of genetic programming, developed a new parameter (i.e., the G p criterion) to measure the glass-forming ability of BMGs. It has been suggested that the stability of the liquid phase and the glass's resistance to crystallization should be considered simultaneously to increase the GFA of alloys 17 . Therefore, it is impossible to use some of these parameters as indicators of GFA. For example, the ΔT x parameter only considers the stability of the glassy phase and not the ease of glass formation. Similarly, T rg (= T g /T l ) is a parameter that only considers the ease of glass formation 18 . In addition, although developing different parameters able scientists to evaluate and predict the glass-forming ability of alloy systems, most of these expressions have a low correlation with the glass-forming ability and will fail in predicting the D max of newly developed alloy systems 16,19 .
In recent years, ML has been known as one of the best routes in predicting the glass-forming ability of BMGs [20][21][22][23] . Scientists have used ML as a promising route for materials design and discovery 24 . The ML is the science of developing models that will become more efficient over time and is going to replace old methods of finding a solution to relate input features to output features with new attitudes to solving the problem with artificial intelligence 25 . By using ML modeling, scientists are enable to relate characteristic features of material

Methodology
The whole process of modeling consists of four steps. The schematic demonstration of these steps is represented in Fig. 1. The very first step of modeling is collecting a suitable dataset. The next step is processing the dataset with input features to attain the desired output. It then can be followed by representing results. The final step is the validation of the model. Data collection. In the ML, it is essential to use a suitable dataset. Due to the purpose of ML and attaining a reliable result from modeling, the dataset used in this article was collected from 44 different published articles 12,19, , including 715 different BMGs. Furthermore, 17 different parameters were used as input features to predict the maximum critical diameter (D max ) of BMGs. It is noticeable that all parameters can be calculated using characteristic temperatures. Although liquid fragility is another important parameter that has a straightforward relationship with the GFA, researchers need to have the relaxation time of metallic glasses to measure this parameter. In most papers used in the dataset, there is insufficient information for determining this parameter. Besides that, since one of the parameters we used as an input feature (i.e., the reduced glass-transition temperature or T rg ) is essential in determining the liquid fragility parameter, adding the liquid fragility parameter will not benefit the model's accuracy 72 . Therefore, the liquid fragility parameter has not been used in this   where the X new is the normalized value of x and max(x) and min(x) are maximum and minimum values of feature x, respectively. Recursive Feature Elimination (RFE) is used with the Random Forest algorithm for feature selection. The RFE starts with pruning the least essential feature from the set of features and is repeated on the pruned set until the desired number of features to select is eventually reached 88 . Then, three characteristic temperatures were added to the model with the highest correlation value and all features. In the following step, fivefold cross-validation is used to split the dataset into five groups. Four groups form the training set, and one performs as a testing set. The whole process of modeling was repeated 100 times. The mean of all modeling was used to ensure the accuracy and generalization ability of the model 89 . Then, the squared correlation coefficient (R 2 ) and mean absolute error (MAE) were selected to evaluate the accuracy and generalization ability of the RF model. Equations (2) and (3) represent the formula of these methods: where n is the number of samples and f(x i ) and y i represent the predicted and experimental values of the i th sample, respectively.

Validation.
For validation of the model, the characteristic temperatures of the four alloys studied in the previous work 90 were used as inputs. The predicted D max of these alloys was compared with the measured values Table 1. Parameters (GFA criteria) expressed by characterization temperatures used as input features in this article (* means the ideal value has not been measured for this parameter). 11

Results and discussion
Feature selection. The Random Forest Regression (RFR) with 5-fold cross-validation was conducted on the BMGs dataset 100 times to have a reliable model. Then, the Recursive Feature Elimination (RFE) algorithm was used to determine the squared correlation coefficient (R 2 ) of models with the different number of selected features, as shown in Fig. 2. It can be seen from Fig. 2a that the correlation of the best model with a certain number of selected features (i.e., the highest correlation for the model with the specific number of selected features) can be divided into three different parts, where the correlation will maintain the same level at that part. Moreover, it is indicated that models with more than 11 selected features have the highest correlation. In fact, this figure indicates that when we select just one parameter (from 17 parameters) and train the model, the testset accuracy of the best model would be around 88.25%. Similarly, when we choose two parameters and train the model, the accuracy will increase to approximately 91%. This increase will continue until we select 11 parameters. At this step, the accuracy would be around 92.5%. After that, the change in accuracy would be imperceptible when we increase the number of selected parameters in training the model. Figure 2b represents the average coefficient of all 100 repeated modeling and illustrates that the model with 13 selected features has the highest average correlation value compared to other models. Although the model with 16 features also had a correlation coefficient near the model with 13 features, it is preferable to have the least number of features 91 . Moreover, overfitting is higher in the models with more parameters 92 . As a result, the model including 13 features was selected as the optimal model for further research in this article.
To ensure that the characteristic temperatures will not be eliminated if we use RFE for all parameters, including characteristic temperatures, we conducted RFE on a model with all 20 parameters. The selected parameters, in this case, were the same as the model with 13 features + 3 characteristic temperatures, i.e., T g , T x, and T l , and the characteristic temperatures were selected as features.
To evaluate the effect of adding characteristic temperatures to the model and thermodynamically guiding the model, T g , T x, and T l were added to the optimal model (i.e., the model with 13 features) and the model with all features. To evaluate every single feature's importance in modeling, the feature importance algorithm is represented in Fig. 3. Figure 3a,b are related to the feature importance of models with 13 different features without and with characteristic temperatures. Also, Fig. 3c,d are for models with all parameters without and with characteristic temperatures. It is evident that in all of the modeling, the parameter of G p is the most crucial feature. There are some reasons why this parameter is so important, and actually, it has the highest importance in all of the modelings. First, the Gp expression has been extracted using genetic programming. Scientists could generate several different solutions using this method and choose the most reliable solution 93 . Secondly, the expression itself can be divided into two different parts that are shown in Eq. (4): Increasing Tg in the first part of Eq. (4) will increase the viscosity, and the glass-formation process will facilitate. Also, T l -T x area is prone to crystallization, so minimizing this part of the Equation is desirable 75 . Besides, it is optimal to have an extended supercooled liquid area (i.e., T x − T g ) in BMGs as crystallization will not occur at this area 78 . On the other hand, the G p expression proposes that high T x will increase liquid stability and lead to better glass-forming ability. Based on Wakasugi et al. 94 , the viscosity of supercooled liquid will increase with (4) www.nature.com/scientificreports/ increasing the ratio of T x /T l, and increasing the viscosity of supercooled liquid will lead to high glass-forming ability. In conclusion, the G p expression inherits the phenomenological attributes of glass-forming ability in BMGs and shows a good correlation with D max 16 .
GFA prediction. Figure 4 shows the testing set's predicted Dmax versus the measured D max . As can be seen, adding characteristic temperatures to the models increases the squared correlation coefficient for the train set and test set and decreases the mean absolute error. The best result is for the model with 13  The effect of characteristic temperatures on GFA prediction. Figure 5 indicates the correlation coefficient for models without characteristic temperatures and those including characteristic temperatures. The R 2 of the test set increased by adding characteristic temperatures in both models, i.e., the model with 13 features and the model with 17 features. It is important to note that mean absolute error was reduced by introducing characteristic temperatures. Characteristic temperatures and selected parameters have a synergic effect on each other because; by using these parameters, the thermodynamic approach of the model will improve. Firstly, all the input parameters are dimensionless relations, including the characteristic temperatures. It is expected that introducing the individual characteristic temperatures will increase the weight of these crucial parameters, which in turn leads the model to be more accurate. In an ideal glassy system with high glass-forming ability, the glass transition temperature (and crystallization temperature) should be very high and the liquidus temperature very low 3 . Secondly, the two critical criteria for high glass-forming ability are the liquid phase's stability and crystallization resistance 76 . Although in most of the parameters, these two factors have been considered, some of them, such as T rg and ΔT x, are only based on the stability of the liquid phase 15,73 . It has been indicated that the stability in the equilibrium state and in the undercooled state depends on the liquidus and glass transition temperatures, respectively. On the other hand, resistance to crystallization is a parameter that has a linear relationship with the glass transition temperature 13,76 . As a result, by introducing the characteristic temperatures as input parameters, ignorance of parameters like ease of glass formation and resistance to crystallization will be compensated.   www.nature.com/scientificreports/ Overfitting is one of the most critical problems in ML modeling that cause us to have a low coefficient correlation in the testing set while having a very well-fitted model on the training set. The Overfitting/Underfitting can be determined using the squared correlation coefficient of the training and test sets 97 . For evaluating the effect of the introduction of characteristic temperatures on the overfitting problem 24 , the number of models with a difference of R 2 value between the training set and testing set in the range of 1-5% and 5-10% are plotted respectively as shown in Fig. 6. The overfitting problem and this difference are dependent on each other. The overfitting problem is more severe in cases where the difference of R 2 value is higher 98 . It is evident that adding the characteristic temperatures to the models reduces the occurrence of overfitting and improves the efficiency and accuracy of the model.

Validation.
In order to validate the established model (the one including 13 features and characteristic temperatures), the measured D max as the glass-forming ability of four different alloys along with the related data were extracted from reference 90 . The data of alloys entered into the model, and then the predicted D max is compared with the experimental results, as shown in Table 2. The modeling results were very close to the experimental results, as shown by low discrepancy between the experimental and predicted results, which is in the range of regression value (Fig. 4b). This works shows that neither 13 dimensionless parameters proposed by different researchers nor three characteristic temperatures proposed in this work could predict the GFA accurately. However, when three characteristic temperatures plus 13 parameters were employed by ML modeling, a rather more accurate result was depicted. In fact, most of the GFA results in literature have been obtained by trial-and-error experiments without sufficient scientific background, e.g., thermodynamic analysis. This matter is mostly due to non-equilibrium and complicated feathers of solidification in BMG alloys. The successful introduction of three characteristic temperatures, which are closely related to the thermodynamics of BMG alloys, bears in mind that more scientific exploring on the composition-structure-properties of BMG could be very informative to provide enough knowledge for predicting GFA to lessen too many experimental studies.

Conclusion
In summary, this work conducted a random forest model combined with Recursive Feature Elimination to predict the D max (GFA) of BMGs on a dataset of 714 different alloys. The results showed that the model with 13 selected features has the highest coefficient correlation value (R 2 = 92%) and is very efficient. The feature importance algorithm showed that the G p parameter is the most crucial parameter in all modeling. Adding characteristic temperatures, i.e., T g , T x, and T l , to the models will improve the accuracy (R 2 up to 95%) and efficiency of models, resulting in overcoming the overfitting problem. Comparing the measured and predicted values of D max for a set of newly developed BMGs indicated that the model is reliable and can be used for predicting the GFA of BMGs.

Data availability
The collected dataset is represented in "Supplementary Information Files". Developed codes for Machine-Learning are accessible through http:// github. com/ naha7 789/ GFA-predi cting.